%% master_baseline_CES_approx.m

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% AUTHOR: Produced for by JDG/IM based on initial code for cobb-douglas
% production in "Firm Model" directory.
%
% DATE: January 2019.
%
% Purpose: This is the main file that is used to execute all of the Matlab
% files for the initial convergence loop that will provide you with a 
% Matlab MAT file that can then be used for the counterfactual exercises. 
% The convergence loop here is found by setting
% all shocks to 1 (no shocks) and the wedges to zero (looping to a 'new'
% fixed point).
%
%  NEW: crucial ingredient is nested-CES production functions. This implies
%  both intermediate input and labor shares have to be updated in the hat
%  algebra, wherase the baseline treated these as data. This in turn
%  implies four variables to be updated in t+1 (see page 32 of notes):
%   (1) total cost function, b
%   (2) price of intermediate bundle, P^m
%   (3) labor share, pi^l
%   (4) intermiedate input share, pi^M
%
% Four files are run:
%  1) Step_1_Import_WIOD.m: import WIOD sector-level data.
%  2) Step_2_Real_Firms.m: import French firm-level data.
%  3) Step_3_Market_Clearing_Expenditure_CES.m: solve for "wedges" that
%     allow model to match observed trade data in levels (steady-state).
%  4) Step_4_Algorithm_baseline_CES.m: solve for changes (hat algebra) in 
%     order to eliminate wedges. Then use changes to generate adjusted
%     steady-state level variables to perform counterfactuals on. 
%
% Output: MAT/Step_4_setup_`alphaT''rhoT'.mat
% 
% Note: "T" refers to type, and appended given based on parameter names we
%       use for values in the coding. Not all will be used depending on
%       type of file being generated in order to keep names shorter.
%
% alphaT: setting of labor share from data (e.g., homog or hetero)
%
% rhoT: elasticity of substitution in deamdn for goods within a sector
%
% lambdaT: elasticity of subtitution between labor and intermediate bundle
%
% etaT: elastictiy of substitution between intermeidate goods in bundle
%
% shockT: type of shock scheme applied       
%   
% Finally: for all code, we will change the subscripts to match those in
% master_notes3.pdf and the draft. I.e., we will change to {mn,ij} pairs,
% which represent country pair {mn} and sector pair {ij}.
%
% Last Updated: 24/01/2019 by JDG
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



global alphaT rhoT lambdaT etaT rho sigma sigmaT lambda eta varphi ...
    tol maxit Pi_hat_n D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    N J FI_J france countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO ... 
    sL_n sPi_n sD_n alpha_WIOD_j alpha_WIOD_francej labour_share_PWT
% shockT oneminus_alpha_vector  sector_algorithm
    

%% SET PARAMETERS VALUES %%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% alphaT: 
%
% 1 = PWT labor share for all countries and Labor share for French firms 
%     are at the COUNTRY level.
%
% 2 = PWT labor share for all countries and WIOD scaled labor share for 
%     French firms. Labor share for French firms are at the FIRM level.
%
% 3 = WIOD labor share for all SECTORS with the labor share at the SECTOR 
%     level for French firms.
%
% 4 = WIOD labor share for all SECTORS with the labor share at the FIRM 
%     level for French firms 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%alphaT = 4;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rhoT:
%
% Value for rho. Will apply to J number of sectors to create a vector of
% homogeneous rhos across sectors.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%rhoT = '3';

% Number of sectors we use:

J = 32;

%rho = 3*ones(J,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Value for sigma
%
% This is the elasticity of substitution across sector consumption. Given
% we do not adjust this anymore after running C-D production we set to the
% value we think is most reasonable: 1.5. There is therefore no longer a
% "type" variable.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% sigma = 3*ones(J,1);
% sigmaT = '3';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% lambdaT:
%
% Value for lambda. We might experiment on this so have a type.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%lambdaT = '1.000001';

%lambda = 1.000001;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% etaT:
%
% Value for eta. We might experiment on this so have a type.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%etaT = '1.000001';
%eta = 1.000001;

% Frisch elasticity for GHH preferences

%varphi = 3;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% shockT: not needed for baseline
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% RUNNING THE ALGORITHM %%

disp('BEGINNING SETUP')

% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Step_1_Import_WIOD95:    
% %
% % Imports the WIOD data, aggregates countries and removes sectors that we 
% % will not be analysing. 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
 Step_1_Import_WIOD_CES
% 
 save MAT/WIOD.mat france N J gamma start PC_mnj PC_mj PC_nj PC_n ...
            pi_c_mnj pi_c_nj Z_mnj X_mnj X_nj pi_mnj labour_share_PWT...
            countrysecD X_mn sum_by_country_dummy alpha_WIOD_j ...
            alpha_WIOD_francej alpha_WIOD_nj;
% 
 disp('STEP 1 COMPLETED SUCCESFULLY')
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Step_2_Real_Firms: 
% % Imports the raw French firm gammas and pi. Then reorders the firm data so
% % that it matches the WIOD ordering and save this data for use later
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
 Step_2_Real_Firms_CES
% 
 Step_2_save =[' GAMMA_sorted_matrix alpha_f_sorted pi_mnj_sorted '...
      ' firmsecD_sorted ISO start_sorted FI_J firms_per_sector_sorted'...
      ' VA_f_Stata turnover_f_Stata turnover_f_Stata_share WIOD_sectors' ...
      ' sector SECTORS_sorted '];
% 
 eval([strcat('save MAT/Step_2_firm_data_alpha',num2str(alphaT),...
      '.mat',Step_2_save, ' -v7.3')]);
% 
 disp('STEP 2 COMPLETED SUCCESFULLY')
% 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Step_3_Market_Clearing_Expenditure: 
% %
% % Finds the market clearing wedges. These are the wedges needed to have 
% % the market clear. 
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
 Step_3_Market_Clearing_Expenditure_CES
% 
%  % Assigning the variables we wish to save
% 
% 
 Step_3_save= [' zeta_mnj'];
% 
%  % Saving the output
% 
 eval([strcat('save MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_3_alpha',num2str(alphaT),...
      '_rho',rhoT,'.mat',Step_3_save,' -v7.3;')]);
% 
 disp('STEP 3 COMPLETED SUCCESFULLY')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Step_4_Algorithm: 
% 
% Implements the convergence algorthim using hat algebra. Such that a 
% X_hat = X_{t+1}/X_t. Here we calculate the new  hat variables 
% (pseudo changes) for when there are no wedges in the model.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Load in the data

load MAT/WIOD.mat
eval([strcat('load MAT/Step_2_firm_data_alpha',num2str(alphaT),'.mat')]);
eval([strcat('load MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_3_alpha',num2str(alphaT),...
    '_rho',rhoT,'.mat')]);

Step_4_Algorithm_baseline_CES_approx

% Updating the data based on the results of the baseline algorithm

pi_M = pi_M1;
pi_M_f = pi_M_f1;
pi_l = pi_l1;
pi_l_f = pi_l_f1;
pi_mnj = pi_mnj1;
pi_mnjf = pi_mnjf1;
pi_c_mnj = pi_c_mnj1;
pi_c_nj = pi_c_nj1;
zeta_mnj = zeros(J*N,N);
X_mnj = X_hat_mnj.*X_mnj;
PC_n = PC_hat_n1.*PC_n;  
final_goods = pseudo_final_goods1;
intermediates = pseudo_intermediates1;

% Updating sectoral pi_M and pi_l matrix for France so that they match with
% the firm-level ones
% Sectoral alphas are a weighted average of firm-level ones based on sales
% weights
X_nkif = firmsecD_sorted*X_mnj(start(france):start(france+1)-1,:);
pi_X_nkif=pi_mnjf.*X_nkif;
pi_X_nif=sum(pi_X_nkif,2);
pi_X_ni=firmsecD_sorted'*pi_X_nif;
pi_X_ni_ = firmsecD_sorted*pi_X_ni;
pi_l_pi_X_nif=pi_l_f.*pi_X_nif./pi_X_ni_;
pi_l_sect = firmsecD_sorted'*pi_l_pi_X_nif;
pi_l_sect_old = pi_l(start(france):start(france+1)-1,:);
corr(pi_l_sect,pi_l_sect_old)
pi_l(start(france):start(france+1)-1,:) = pi_l_sect;

% Sectoral gammas are a weighted average of firm-level ones based on input
% purchases weights
one_minus_alpha = 1-pi_l_f;
one_minus_alpha_pi_X_nif = one_minus_alpha.*pi_X_nif;
one_minus_alpha_pi_X_ni=firmsecD_sorted'*one_minus_alpha_pi_X_nif;
one_minus_alpha_pi_X_ni_ = firmsecD_sorted*one_minus_alpha_pi_X_ni;
one_minus_alpha_pi_X_ni_repmat = repmat(one_minus_alpha_pi_X_ni_,[1,N*J]);
one_minus_alpha_pi_X_nif_repmat = repmat(one_minus_alpha_pi_X_nif,[1,N*J]);
pi_M_f_weight_I=pi_M_f.*one_minus_alpha_pi_X_nif_repmat./one_minus_alpha_pi_X_ni_repmat;
pi_M_sect = firmsecD_sorted'*pi_M_f_weight_I;
pi_M_sect_old = pi_M(:,start(france):start(france+1)-1)';
pi_M_sect_reshape=reshape(pi_M_sect,[],1);
pi_M_sect_old_reshape=reshape(pi_M_sect_old,[],1);
corr(pi_M_sect_reshape,pi_M_sect_old_reshape)
corr(pi_M_sect_reshape(start(france):start(france+1)-1,:),pi_M_sect_old_reshape(start(france):start(france+1)-1,:))
pi_M(:,start(france):start(france+1)-1) = pi_M_sect';

% Saving the baseline calibration, the fixed point is the starting point
% for the counterfactual exercises. This needs to be done as we need to run
% a setup initially where there are no wedges (market clearing condition
% holds), this calculates the new X_t+1 and the other level variables.

Step_4_save = [' pi_M pi_M_f pi_l pi_l_f pi_mnj pi_mnjf pi_c_mnj pi_c_nj ' ...
               ' zeta_mnj X_mnj PC_n final_goods intermediates'];

 eval([strcat('save MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_setup_alpha',num2str(alphaT),...
     '_rho',rhoT,'_approx.mat',Step_4_save,' -v7.3;')]);

disp('STEP 4 BASELINE COMPLETED SUCCESFULLY')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


disp('SETUP COMPLETED SUCCESFULLY')